
#ifndef __projection_errors_h__
#define __projection_errors_h__

#include "petsc.h"
#include "libmesh.h"
#include "mesh.h"
#include "equation_systems.h"
#include "petsc_vector.h"
#include "perf_log.h"

#include "point.h"
#include "function_base.h"

#define PRIMARY_SYSTEM_NAME "Implicit_System"
#define AUXILLARY_SYSTEM_NAME "Explicit_System"

// declare our functions
void generate_mesh(unsigned int dim, unsigned int ps, const std::string& order, Mesh& mesh) ;

//int libmesh_project_solution(EquationSystems& src, EquationSystems& target, unsigned int dim=3, bool verbose=false, std::string output_file="");

//int moab_project_solution(EquationSystems& src, EquationSystems& target, unsigned int dim=3, bool verbose=false, std::string output_file="");

void fill_auxillary_systems(EquationSystems& es, std::string aux_sys_name) ;

// Exact solution function prototype.
Real exact_solution (const Real x,
                     const Real y = 0.,
                     const Real z = 0.);

void compute_errors(EquationSystems& es, std::vector<Real>& errors, std::string primary_sys_name) ;

/**
 * This is the exact solution functor.
 */
class ExactPoisson : public FunctionBase<Real>
{
public:
  ~ExactPoisson() { }

  virtual Real operator() (const Point &p, const Real =0.)
  {
    return exact_solution (p(0), p(1), p(2)) ;
  }

  virtual void operator() (const Point &p, const Real , DenseVector<Real> & 	output)
  {
    Real exactval = exact_solution (p(0), p(1), p(2)) ;
    for (unsigned int i=0; i < output.size(); ++i)
      output(i) = exactval;
  }

  AutoPtr< FunctionBase<Real> > clone	()	 const
  {
    return AutoPtr<FunctionBase<Real> > (new ExactPoisson()) ;
  }

};

class ProjectionBase
{
public:
  ProjectionBase(bool is_verbose) : verbose(is_verbose)
  {
    projection_log = NULL;
    source = NULL;
    target = NULL;
  }

  virtual ~ProjectionBase()
  {
    if (projection_log)
      delete projection_log ;
    source = NULL;
    target = NULL;
  }

  void set_source_equation_system(EquationSystems& es) ;

  void set_target_equation_system(EquationSystems& es) ;

  virtual std::string type() = 0 ;

  virtual void init() ;

  virtual int apply(std::string output_file="") = 0 ;

  PerfLog* projection_log ;

protected:

  EquationSystems* source ;
  EquationSystems* target ;
  const bool verbose;

};

// inline methods

inline
void ProjectionBase::set_source_equation_system(EquationSystems& es)
{
  source = &es;
}


inline
void ProjectionBase::set_target_equation_system(EquationSystems& es)
{
  target = &es;
}


inline
void ProjectionBase::init()
{
  libmesh_assert(this->source != NULL) ;
  libmesh_assert(this->target != NULL) ;

  this->projection_log = new PerfLog (this->type());
}



#endif

